/*==============================================================================
SIPP1991数据集DDML分析 - 增强版（含图表绘制和Word导出）
================================================================================

研究目标：
使用双重机器学习(Double/Debiased Machine Learning, DDML)方法分析401(k)参与
对净金融资产的影响，控制混杂因素后获得因果效应估计。

数据来源：
1991年收入与计划参与调查(Survey of Income and Program Participation, SIPP)

变量说明：
- net_tfa: 净金融资产（千美元）- 结果变量Y
- e401: 401(k)计划参与（1=参与，0=未参与）- 处理变量D  
- age: 年龄
- inc: 收入（千美元）
- educ: 教育年限
- fsize: 家庭规模
- marr: 婚姻状况（1=已婚，0=未婚）
- twoearn: 双收入家庭（1=是，0=否）
- db: 固定收益计划参与（1=是，0=否）
- pira: 个人退休账户参与（1=是，0=否）
- hown: 住房所有权（1=是，0=否）

分析方法：
1. 描述性统计分析
2. 数据可视化
3. DDML因果效应估计
4. 结果解释和报告生成

作者：Stata ML Course
日期：2025-11-06
==============================================================================*/

clear all
set more off
capture log close

* 创建输出目录
capture mkdir "output/sipp1991"
capture mkdir "output/sipp1991/figures"
capture mkdir "output/sipp1991/tables"

* 开始日志记录
log using "output/sipp1991/sipp1991_ddml_analysis.log", replace

/*------------------------------------------------------------------------------
第一部分：数据加载和描述性统计
------------------------------------------------------------------------------*/

* 加载数据
use "data/sipp1991.dta", clear

* 设置全局变量
global Y net_tfa
global D e401
global X age inc educ fsize marr twoearn db pira hown

* 数据清理
foreach var of varlist $Y $D $X {
    drop if missing(`var')
}

* 描述性统计
summarize $Y $D $X, detail

* 按401(k)参与状态分组统计
tabstat $Y, by($D) stats(mean sd n)

/*------------------------------------------------------------------------------
第二部分：数据可视化
------------------------------------------------------------------------------*/

* 2.1 净金融资产分布直方图
histogram $Y, width(5) frequency ///
    title("净金融资产分布") ///
    xtitle("净金融资产（千美元）") ytitle("频数") ///
    note("数据来源：SIPP 1991") scheme(s2color)
graph export "output/sipp1991/figures/net_tfa_histogram.png", as(png) replace

* 2.2 按401(k)参与状态的净金融资产对比
graph bar (mean) $Y, over($D) ///
    title("按401(k)参与状态的平均净金融资产") ///
    ytitle("平均净金融资产（千美元）") blabel(bar, format(%9.1f)) ///
    legend(label(1 "未参与") label(2 "参与")) scheme(s2color)
graph export "output/sipp1991/figures/net_tfa_by_e401.png", as(png) replace

* 2.3 净金融资产箱线图
graph box $Y, over($D) ///
    title("净金融资产分布箱线图") ytitle("净金融资产（千美元）") ///
    legend(label(1 "未参与") label(2 "参与")) scheme(s2color)
graph export "output/sipp1991/figures/net_tfa_boxplot.png", as(png) replace

* 2.4 年龄与净金融资产关系
twoway (scatter $Y age, mcolor(blue%40)) ///
       (lfit $Y age, lcolor(red) lwidth(thick)), ///
    title("年龄与净金融资产关系") ///
    xtitle("年龄") ytitle("净金融资产（千美元）") ///
    legend(order(1 "数据点" 2 "线性拟合")) scheme(s2color)
graph export "output/sipp1991/figures/net_tfa_age_scatter.png", as(png) replace

* 2.5 收入与净金融资产关系
twoway (scatter $Y inc, mcolor(green%40)) ///
       (lfit $Y inc, lcolor(red) lwidth(thick)), ///
    title("收入与净金融资产关系") ///
    xtitle("收入（千美元）") ytitle("净金融资产（千美元）") ///
    legend(order(1 "数据点" 2 "线性拟合")) scheme(s2color)
graph export "output/sipp1991/figures/net_tfa_inc_scatter.png", as(png) replace

* 2.6 教育与净金融资产关系
graph bar (mean) $Y, over(educ) ///
    title("不同教育水平的平均净金融资产") ///
    ytitle("平均净金融资产（千美元）") blabel(bar, format(%9.1f)) scheme(s2color)
graph export "output/sipp1991/figures/net_tfa_education.png", as(png) replace

* 2.7 婚姻状况与净金融资产
graph bar (mean) $Y, over(marr) ///
    title("婚姻状况与净金融资产") ///
    ytitle("平均净金融资产（千美元）") blabel(bar, format(%9.1f)) ///
    legend(label(1 "未婚") label(2 "已婚")) scheme(s2color)
graph export "output/sipp1991/figures/net_tfa_marriage.png", as(png) replace

* 2.8 住房所有权与净金融资产
graph bar (mean) $Y, over(hown) ///
    title("住房所有权与净金融资产") ///
    ytitle("平均净金融资产（千美元）") blabel(bar, format(%9.1f)) ///
    legend(label(1 "无房") label(2 "有房")) scheme(s2color)
graph export "output/sipp1991/figures/net_tfa_homeownership.png", as(png) replace

/*------------------------------------------------------------------------------
第三部分：DDML分析
------------------------------------------------------------------------------*/

* 设置随机种子
set seed 42

* 初始化DDML
ddml init partial, kfolds(5) reps(5)

* 估计结果变量Y的条件期望 - 使用多种方法
ddml E[Y|X]: regress $Y $X
ddml E[Y|X]: pystacked $Y $X, type(reg) method(rf)
ddml E[Y|X]: pystacked $Y $X, type(reg) method(lasso)

* 估计处理变量D的条件期望 - 使用多种方法
ddml E[D|X]: logit $D $X
ddml E[D|X]: pystacked $D $X, type(reg) method(rf)
ddml E[D|X]: pystacked $D $X, type(reg) method(lasso)

* 显示模型描述
ddml desc

* 交叉拟合
ddml crossfit

* 估计因果效应
ddml estimate, robust

* 保存DDML结果
estimates store ddml_results

* 提取DDML估计结果
matrix ddml_b = e(b)
matrix ddml_V = e(V)
local ddml_coef = ddml_b[1,1]
local ddml_se = sqrt(ddml_V[1,1])
local ddml_t = `ddml_coef' / `ddml_se'
local ddml_p = 2*ttail(e(N)-1, abs(`ddml_t'))
local ddml_ci_lower = `ddml_coef' - 1.96*`ddml_se'
local ddml_ci_upper = `ddml_coef' + 1.96*`ddml_se'

* 运行OLS回归作为对比
quietly regress $Y $D $X
estimates store ols_results
local ols_coef = _b[$D]
local ols_se = _se[$D]
local ols_t = `ols_coef' / `ols_se'
local ols_p = 2*ttail(e(df_r), abs(`ols_t'))
local ols_ci_lower = `ols_coef' - 1.96*`ols_se'
local ols_ci_upper = `ols_coef' + 1.96*`ols_se'
local ols_n = e(N)
local ols_r2 = e(r2)

* 运行所有组合的DDML估计
ddml estimate, robust allcombos

/*------------------------------------------------------------------------------
第四部分：结果可视化
------------------------------------------------------------------------------*/

* 4.1 创建效应估计图
preserve
clear
input str20 method estimate lower upper
"DDML估计" `ddml_coef' `ddml_ci_lower' `ddml_ci_upper'
"OLS估计" `ols_coef' `ols_ci_lower' `ols_ci_upper'
end

gen method_num = _n
twoway (rcap lower upper method_num, horizontal lcolor(blue)) ///
       (scatter method_num estimate, mcolor(red) msize(large)), ///
    title("401(k)参与对净金融资产的因果效应估计") ///
    subtitle("正值表示参与增加净金融资产") ///
    xtitle("效应估计值（千美元）") ytitle("估计方法") ///
    ylabel(1 "DDML" 2 "OLS", valuelabel) ///
    xline(0, lpattern(dash) lcolor(gray)) legend(off) scheme(s2color)
graph export "output/sipp1991/figures/causal_effect_comparison.png", as(png) replace
restore

* 4.2 创建预测值vs实际值图（如果可用）
capture predict yhat, xb
if _rc == 0 {
    twoway (scatter $Y yhat, mcolor(blue%40)) ///
           (function y=x, range(0 100) lcolor(red) lpattern(dash)), ///
        title("预测值vs实际值") ///
        xtitle("预测净金融资产（千美元）") ytitle("实际净金融资产（千美元）") ///
        legend(order(1 "数据点" 2 "完美预测线")) scheme(s2color)
    graph export "output/sipp1991/figures/prediction_vs_actual.png", as(png) replace
}

* 4.3 创建残差图
capture predict residuals, residuals
if _rc == 0 {
    twoway (scatter residuals yhat, mcolor(blue%40)) ///
           (yline(0, lcolor(red) lpattern(dash))), ///
        title("残差图") ///
        xtitle("预测值") ytitle("残差") ///
        legend(off) scheme(s2color)
    graph export "output/sipp1991/figures/residual_plot.png", as(png) replace
}

/*------------------------------------------------------------------------------
第五部分：生成Word报告
------------------------------------------------------------------------------*/

* 保存样本量到宏
quietly count
local total_obs = r(N)

* 创建Word文档
putdocx begin, header(main_header) footer(main_footer)

* 设置页眉
putdocx paragraph, toheader(main_header) font("微软雅黑", 10)
putdocx text ("SIPP1991数据集DDML分析报告 | Double Machine Learning Analysis")

* 设置页脚
putdocx paragraph, tofooter(main_footer) halign(center) font("微软雅黑", 9)
putdocx text ("第 ")
putdocx pagenumber
putdocx text (" 页")

* 标题页
putdocx paragraph, style(Title) halign(center)
putdocx text ("401(k)参与对净金融资产的影响"), font("微软雅黑", 28, black) bold

putdocx paragraph, halign(center) spacing(after, 10)
putdocx text ("基于双重机器学习的因果效应分析"), font("微软雅黑", 16, "gray") italic

putdocx paragraph, halign(center) spacing(after, 20)
putdocx text ("生成日期: "), font("微软雅黑", 11)
putdocx text ("`c(current_date)'"), font("微软雅黑", 11) bold

putdocx paragraph, halign(center)
putdocx text ("数据来源：SIPP 1991"), font("微软雅黑", 11)

putdocx paragraph, halign(center) spacing(after, 30)
putdocx text ("样本量: `total_obs' 个观测"), font("微软雅黑", 11)

putdocx pagebreak

* 执行摘要
putdocx paragraph, style(Heading1)
putdocx text ("执行摘要"), font("微软雅黑", 18) bold

putdocx textblock begin
本报告使用双重机器学习(DDML)方法分析了401(k)退休计划参与对个人净金融资产的因果效应。
通过控制年龄、收入、教育、家庭规模、婚姻状况、双收入状况、固定收益计划参与、
个人退休账户参与和住房所有权等混杂因素，我们获得了401(k)参与对净金融资产影响的无偏估计。
putdocx textblock end

putdocx pagebreak

* 第一章：数据概览
putdocx paragraph, style(Heading1)
putdocx text ("一、数据概览"), font("微软雅黑", 18) bold

putdocx paragraph, style(Heading2) spacing(before, 10)
putdocx text ("1.1 数据集基本信息"), font("微软雅黑", 14) bold

putdocx textblock begin
本研究使用1991年收入与计划参与调查(SIPP)数据集，包含`total_obs'个个人观测值。
主要变量包括净金融资产、401(k)参与状态以及一系列控制变量。
putdocx textblock end

* 数据基本信息表
putdocx table info = (6, 2), border(all, single, black) layout(autofitcontents)
putdocx table info(1,1) = ("项目"), bold font("微软雅黑", 11)
putdocx table info(1,2) = ("数值"), bold font("微软雅黑", 11)
putdocx table info(1,.), shading("lightblue")

putdocx table info(2,1) = ("总观测数"), font("微软雅黑", 10)
putdocx table info(2,2) = ("`total_obs'"), font("微软雅黑", 10)

putdocx table info(3,1) = ("结果变量"), font("微软雅黑", 10)
putdocx table info(3,2) = ("净金融资产（千美元）"), font("微软雅黑", 10)

putdocx table info(4,1) = ("处理变量"), font("微软雅黑", 10)
putdocx table info(4,2) = ("401(k)参与（0/1）"), font("微软雅黑", 10)

putdocx table info(5,1) = ("控制变量数"), font("微软雅黑", 10)
putdocx table info(5,2) = ("9个"), font("微软雅黑", 10)

putdocx table info(6,1) = ("分析方法"), font("微软雅黑", 10)
putdocx table info(6,2) = ("DDML"), font("微软雅黑", 10)

putdocx pagebreak

* 第二章：描述性统计分析
putdocx paragraph, style(Heading1)
putdocx text ("二、描述性统计分析"), font("微软雅黑", 18) bold

putdocx paragraph, style(Heading2) spacing(before, 10)
putdocx text ("2.1 全样本描述性统计"), font("微软雅黑", 14) bold

putdocx textblock begin
表1展示了所有主要变量的描述性统计信息，包括观测数、均值、标准差、最小值和最大值。
putdocx textblock end

* 表1: 全样本描述性统计（三线表）
quietly {
    * 收集所有变量的统计信息
    summarize net_tfa
    local nfa_n = r(N)
    local nfa_mean = r(mean)
    local nfa_sd = r(sd)
    local nfa_min = r(min)
    local nfa_max = r(max)

    summarize e401
    local e401_n = r(N)
    local e401_mean = r(mean)
    local e401_sd = r(sd)
    local e401_min = r(min)
    local e401_max = r(max)

    summarize age
    local age_n = r(N)
    local age_mean = r(mean)
    local age_sd = r(sd)
    local age_min = r(min)
    local age_max = r(max)

    summarize inc
    local inc_n = r(N)
    local inc_mean = r(mean)
    local inc_sd = r(sd)
    local inc_min = r(min)
    local inc_max = r(max)

    summarize educ
    local educ_n = r(N)
    local educ_mean = r(mean)
    local educ_sd = r(sd)
    local educ_min = r(min)
    local educ_max = r(max)

    summarize fsize
    local fsize_n = r(N)
    local fsize_mean = r(mean)
    local fsize_sd = r(sd)
    local fsize_min = r(min)
    local fsize_max = r(max)
}

putdocx paragraph
putdocx text ("表1 全样本描述性统计"), font("微软雅黑", 10) bold

* 创建三线表
putdocx table tbl1 = (8, 6), border(all, nil) layout(autofitcontents)
* 顶部粗线
putdocx table tbl1(1,.), border(top, single, black, 1.5)
putdocx table tbl1(2,.), border(bottom, single, black, 1)
* 底部粗线
putdocx table tbl1(8,.), border(bottom, single, black, 1.5)

* 表头
putdocx table tbl1(1,1) = ("变量"), font("微软雅黑", 10) bold halign(left)
putdocx table tbl1(1,2) = ("观测数"), font("微软雅黑", 10) bold halign(center)
putdocx table tbl1(1,3) = ("均值"), font("微软雅黑", 10) bold halign(center)
putdocx table tbl1(1,4) = ("标准差"), font("微软雅黑", 10) bold halign(center)
putdocx table tbl1(1,5) = ("最小值"), font("微软雅黑", 10) bold halign(center)
putdocx table tbl1(1,6) = ("最大值"), font("微软雅黑", 10) bold halign(center)

* 数据行
putdocx table tbl1(2,1) = ("净金融资产（千美元）"), font("微软雅黑", 9) halign(left)
putdocx table tbl1(2,2) = ("`nfa_n'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc)
putdocx table tbl1(2,3) = ("`nfa_mean'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(2,4) = ("`nfa_sd'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(2,5) = ("`nfa_min'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(2,6) = ("`nfa_max'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)

putdocx table tbl1(3,1) = ("401(k)参与（0/1）"), font("微软雅黑", 9) halign(left)
putdocx table tbl1(3,2) = ("`e401_n'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc)
putdocx table tbl1(3,3) = ("`e401_mean'"), font("微软雅黑", 9) halign(center) nformat(%9.3f)
putdocx table tbl1(3,4) = ("`e401_sd'"), font("微软雅黑", 9) halign(center) nformat(%9.3f)
putdocx table tbl1(3,5) = ("`e401_min'"), font("微软雅黑", 9) halign(center) nformat(%9.0f)
putdocx table tbl1(3,6) = ("`e401_max'"), font("微软雅黑", 9) halign(center) nformat(%9.0f)

putdocx table tbl1(4,1) = ("年龄（岁）"), font("微软雅黑", 9) halign(left)
putdocx table tbl1(4,2) = ("`age_n'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc)
putdocx table tbl1(4,3) = ("`age_mean'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(4,4) = ("`age_sd'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(4,5) = ("`age_min'"), font("微软雅黑", 9) halign(center) nformat(%9.0f)
putdocx table tbl1(4,6) = ("`age_max'"), font("微软雅黑", 9) halign(center) nformat(%9.0f)

putdocx table tbl1(5,1) = ("收入（千美元）"), font("微软雅黑", 9) halign(left)
putdocx table tbl1(5,2) = ("`inc_n'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc)
putdocx table tbl1(5,3) = ("`inc_mean'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(5,4) = ("`inc_sd'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(5,5) = ("`inc_min'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(5,6) = ("`inc_max'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)

putdocx table tbl1(6,1) = ("教育年限（年）"), font("微软雅黑", 9) halign(left)
putdocx table tbl1(6,2) = ("`educ_n'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc)
putdocx table tbl1(6,3) = ("`educ_mean'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(6,4) = ("`educ_sd'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(6,5) = ("`educ_min'"), font("微软雅黑", 9) halign(center) nformat(%9.0f)
putdocx table tbl1(6,6) = ("`educ_max'"), font("微软雅黑", 9) halign(center) nformat(%9.0f)

putdocx table tbl1(7,1) = ("家庭规模"), font("微软雅黑", 9) halign(left)
putdocx table tbl1(7,2) = ("`fsize_n'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc)
putdocx table tbl1(7,3) = ("`fsize_mean'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(7,4) = ("`fsize_sd'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl1(7,5) = ("`fsize_min'"), font("微软雅黑", 9) halign(center) nformat(%9.0f)
putdocx table tbl1(7,6) = ("`fsize_max'"), font("微软雅黑", 9) halign(center) nformat(%9.0f)

* 表格注释
putdocx paragraph
putdocx text ("注：本表展示了主要变量的描述性统计。样本量为`total_obs'个观测。"), font("微软雅黑", 8) italic

putdocx pagebreak

* 2.2 按401(k)参与状态分组统计（三线表）
putdocx paragraph, style(Heading2) spacing(before, 10)
putdocx text ("2.2 按401(k)参与状态分组统计"), font("微软雅黑", 14) bold

putdocx textblock begin
表2展示了按401(k)参与状态分组的净金融资产及控制变量的统计信息，
以便比较参与组和非参与组的差异。
putdocx textblock end

* 计算分组统计
quietly {
    * 未参与组
    summarize net_tfa if e401 == 0
    local nfa_mean0 = r(mean)
    local nfa_sd0 = r(sd)
    local nfa_n0 = r(N)

    summarize age if e401 == 0
    local age_mean0 = r(mean)
    local age_sd0 = r(sd)

    summarize inc if e401 == 0
    local inc_mean0 = r(mean)
    local inc_sd0 = r(sd)

    summarize educ if e401 == 0
    local educ_mean0 = r(mean)
    local educ_sd0 = r(sd)

    summarize fsize if e401 == 0
    local fsize_mean0 = r(mean)

    summarize marr if e401 == 0
    local marr_mean0 = r(mean)

    summarize twoearn if e401 == 0
    local twoearn_mean0 = r(mean)

    * 参与组
    summarize net_tfa if e401 == 1
    local nfa_mean1 = r(mean)
    local nfa_sd1 = r(sd)
    local nfa_n1 = r(N)

    summarize age if e401 == 1
    local age_mean1 = r(mean)
    local age_sd1 = r(sd)

    summarize inc if e401 == 1
    local inc_mean1 = r(mean)
    local inc_sd1 = r(sd)

    summarize educ if e401 == 1
    local educ_mean1 = r(mean)
    local educ_sd1 = r(sd)

    summarize fsize if e401 == 1
    local fsize_mean1 = r(mean)

    summarize marr if e401 == 1
    local marr_mean1 = r(mean)

    summarize twoearn if e401 == 1
    local twoearn_mean1 = r(mean)

    * 计算差异
    local nfa_diff = `nfa_mean1' - `nfa_mean0'
    local age_diff = `age_mean1' - `age_mean0'
    local inc_diff = `inc_mean1' - `inc_mean0'
    local educ_diff = `educ_mean1' - `educ_mean0'
    local fsize_diff = `fsize_mean1' - `fsize_mean0'
    local marr_diff = `marr_mean1' - `marr_mean0'
    local twoearn_diff = `twoearn_mean1' - `twoearn_mean0'
}

putdocx paragraph
putdocx text ("表2 按401(k)参与状态分组的描述性统计"), font("微软雅黑", 10) bold

* 创建三线表
putdocx table tbl2 = (9, 4), border(all, nil) layout(autofitcontents)
putdocx table tbl2(1,.), border(top, single, black, 1.5)
putdocx table tbl2(2,.), border(bottom, single, black, 1)
putdocx table tbl2(9,.), border(bottom, single, black, 1.5)

* 表头
putdocx table tbl2(1,1) = ("变量"), font("微软雅黑", 10) bold halign(left)
putdocx table tbl2(1,2) = ("未参与组"), font("微软雅黑", 10) bold halign(center)
putdocx table tbl2(1,3) = ("参与组"), font("微软雅黑", 10) bold halign(center)
putdocx table tbl2(1,4) = ("差异"), font("微软雅黑", 10) bold halign(center)

* 数据行
putdocx table tbl2(2,1) = ("净金融资产（千美元）"), font("微软雅黑", 9) halign(left)
putdocx table tbl2(2,2) = ("`nfa_mean0'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl2(2,3) = ("`nfa_mean1'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl2(2,4) = ("`nfa_diff'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)

putdocx table tbl2(3,1) = (""), font("微软雅黑", 8) halign(left)
putdocx table tbl2(3,2) = ("(`nfa_sd0')"), font("微软雅黑", 8) halign(center) italic
putdocx table tbl2(3,3) = ("(`nfa_sd1')"), font("微软雅黑", 8) halign(center) italic
putdocx table tbl2(3,4) = (""), font("微软雅黑", 8) halign(center)

putdocx table tbl2(4,1) = ("年龄（岁）"), font("微软雅黑", 9) halign(left)
putdocx table tbl2(4,2) = ("`age_mean0'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl2(4,3) = ("`age_mean1'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl2(4,4) = ("`age_diff'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)

putdocx table tbl2(5,1) = ("收入（千美元）"), font("微软雅黑", 9) halign(left)
putdocx table tbl2(5,2) = ("`inc_mean0'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl2(5,3) = ("`inc_mean1'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl2(5,4) = ("`inc_diff'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)

putdocx table tbl2(6,1) = ("教育年限（年）"), font("微软雅黑", 9) halign(left)
putdocx table tbl2(6,2) = ("`educ_mean0'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl2(6,3) = ("`educ_mean1'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl2(6,4) = ("`educ_diff'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)

putdocx table tbl2(7,1) = ("已婚比例"), font("微软雅黑", 9) halign(left)
putdocx table tbl2(7,2) = ("`marr_mean0'"), font("微软雅黑", 9) halign(center) nformat(%9.3f)
putdocx table tbl2(7,3) = ("`marr_mean1'"), font("微软雅黑", 9) halign(center) nformat(%9.3f)
putdocx table tbl2(7,4) = ("`marr_diff'"), font("微软雅黑", 9) halign(center) nformat(%9.3f)

putdocx table tbl2(8,1) = ("双收入比例"), font("微软雅黑", 9) halign(left)
putdocx table tbl2(8,2) = ("`twoearn_mean0'"), font("微软雅黑", 9) halign(center) nformat(%9.3f)
putdocx table tbl2(8,3) = ("`twoearn_mean1'"), font("微软雅黑", 9) halign(center) nformat(%9.3f)
putdocx table tbl2(8,4) = ("`twoearn_diff'"), font("微软雅黑", 9) halign(center) nformat(%9.3f)

* 样本量行
putdocx table tbl2(9,1) = ("观测数"), font("微软雅黑", 9) halign(left) bold
putdocx table tbl2(9,2) = ("`nfa_n0'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc) bold
putdocx table tbl2(9,3) = ("`nfa_n1'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc) bold
putdocx table tbl2(9,4) = (""), font("微软雅黑", 9) halign(center)

* 表格注释
putdocx paragraph
putdocx text ("注：括号内为标准差。差异列显示参与组减去未参与组的差值。"), font("微软雅黑", 8) italic

putdocx pagebreak

* 第三章：数据可视化
putdocx paragraph, style(Heading1)
putdocx text ("三、数据可视化分析"), font("微软雅黑", 18) bold

putdocx paragraph, style(Heading2) spacing(before, 10)
putdocx text ("3.1 净金融资产分布"), font("微软雅黑", 14) bold

putdocx textblock begin
下图展示了净金融资产的整体分布情况，可以看出净金融资产呈右偏分布。
putdocx textblock end

putdocx paragraph, halign(center)
putdocx image "output/sipp1991/figures/net_tfa_histogram.png", width(5.5)

putdocx paragraph, style(Heading2) spacing(before, 15)
putdocx text ("3.2 401(k)参与与净金融资产关系"), font("微软雅黑", 14) bold

putdocx textblock begin
下图比较了参与和未参与401(k)计划的个人平均净金融资产。
可以直观看出参与401(k)计划的个人净金融资产明显较高。
putdocx textblock end

putdocx paragraph, halign(center)
putdocx image "output/sipp1991/figures/net_tfa_by_e401.png", width(5.5)

putdocx pagebreak

* 第四章：DDML分析结果
putdocx paragraph, style(Heading1)
putdocx text ("四、DDML因果效应估计"), font("微软雅黑", 18) bold

putdocx paragraph, style(Heading2) spacing(before, 10)
putdocx text ("4.1 模型设定"), font("微软雅黑", 14) bold

putdocx textblock begin
我们使用双重机器学习方法估计401(k)参与对净金融资产的因果效应。
模型设定如下：
- 结果方程：E[Y|X,D] 使用回归、随机森林、Lasso方法
- 处理方程：E[D|X] 使用Logit、随机森林、Lasso方法
- 交叉验证：5折交叉验证，重复5次
- 控制变量：年龄、收入、教育、家庭规模、婚姻状况、双收入状况、固定收益计划、个人退休账户、住房所有权
putdocx textblock end

putdocx paragraph, style(Heading2) spacing(before, 15)
putdocx text ("4.2 因果效应估计结果"), font("微软雅黑", 14) bold

putdocx textblock begin
表3展示了DDML方法和传统OLS方法的估计结果对比。DDML方法通过交叉拟合和去偏处理，
能够更准确地估计因果效应，减少模型选择偏误和过拟合问题。
putdocx textblock end

putdocx paragraph
putdocx text ("表3 401(k)参与对净金融资产的因果效应估计"), font("微软雅黑", 10) bold

* 创建三线表
putdocx table tbl3 = (7, 3), border(all, nil) layout(autofitcontents)
putdocx table tbl3(1,.), border(top, single, black, 1.5)
putdocx table tbl3(2,.), border(bottom, single, black, 1)
putdocx table tbl3(5,.), border(bottom, single, black, 0.5)
putdocx table tbl3(7,.), border(bottom, single, black, 1.5)

* 表头
putdocx table tbl3(1,1) = (""), font("微软雅黑", 10) bold halign(left)
putdocx table tbl3(1,2) = ("DDML估计"), font("微软雅黑", 10) bold halign(center)
putdocx table tbl3(1,3) = ("OLS估计"), font("微软雅黑", 10) bold halign(center)

* 系数估计
putdocx table tbl3(2,1) = ("401(k)效应（千美元）"), font("微软雅黑", 9) halign(left)
putdocx table tbl3(2,2) = ("`ddml_coef'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)
putdocx table tbl3(2,3) = ("`ols_coef'"), font("微软雅黑", 9) halign(center) nformat(%9.2f)

* 标准误
putdocx table tbl3(3,1) = (""), font("微软雅黑", 8) halign(left)
putdocx table tbl3(3,2) = ("(`ddml_se')"), font("微软雅黑", 8) halign(center) italic
putdocx table tbl3(3,3) = ("(`ols_se')"), font("微软雅黑", 8) halign(center) italic

* 95%置信区间
putdocx table tbl3(4,1) = ("95%置信区间"), font("微软雅黑", 9) halign(left)
putdocx table tbl3(4,2) = ("[`ddml_ci_lower', `ddml_ci_upper']"), font("微软雅黑", 8) halign(center)
putdocx table tbl3(4,3) = ("[`ols_ci_lower', `ols_ci_upper']"), font("微软雅黑", 8) halign(center)

* 统计信息
putdocx table tbl3(5,1) = ("观测数"), font("微软雅黑", 9) halign(left)
putdocx table tbl3(5,2) = ("`total_obs'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc)
putdocx table tbl3(5,3) = ("`ols_n'"), font("微软雅黑", 9) halign(center) nformat(%9.0fc)

putdocx table tbl3(6,1) = ("R²"), font("微软雅黑", 9) halign(left)
putdocx table tbl3(6,2) = ("—"), font("微软雅黑", 9) halign(center)
putdocx table tbl3(6,3) = ("`ols_r2'"), font("微软雅黑", 9) halign(center) nformat(%9.3f)

putdocx table tbl3(7,1) = ("控制变量"), font("微软雅黑", 9) halign(left)
putdocx table tbl3(7,2) = ("是"), font("微软雅黑", 9) halign(center)
putdocx table tbl3(7,3) = ("是"), font("微软雅黑", 9) halign(center)

* 表格注释
putdocx paragraph
putdocx text ("注：括号内为稳健标准误。控制变量包括年龄、收入、教育、家庭规模、婚姻状况、双收入状况、固定收益计划、个人退休账户、住房所有权。"), font("微软雅黑", 8) italic
putdocx paragraph
putdocx text ("DDML采用5折交叉验证，重复5次。正值表示参与增加净金融资产。"), font("微软雅黑", 8) italic

putdocx paragraph, style(Heading2) spacing(before, 15)
putdocx text ("4.3 结果可视化"), font("微软雅黑", 14) bold

putdocx textblock begin
图1展示了DDML和OLS两种方法的估计结果及其置信区间的对比。
可以看出，两种方法都显示401(k)参与对净金融资产有显著的正向影响。
putdocx textblock end

putdocx paragraph, halign(center)
putdocx image "output/sipp1991/figures/causal_effect_comparison.png", width(5.5)

putdocx pagebreak

* 4.4 添加完整的回归结果表
putdocx paragraph, style(Heading2) spacing(before, 10)
putdocx text ("4.4 OLS回归详细结果"), font("微软雅黑", 14) bold

putdocx textblock begin
表4展示了OLS回归的完整结果，包括所有控制变量的系数估计。
这有助于理解各个因素对净金融资产的影响。
putdocx textblock end

* 提取OLS所有系数
estimates restore ols_results
local coef_e401 = _b[e401]
local se_e401 = _se[e401]
local coef_age = _b[age]
local se_age = _se[age]
local coef_inc = _b[inc]
local se_inc = _se[inc]
local coef_educ = _b[educ]
local se_educ = _se[educ]
local coef_fsize = _b[fsize]
local se_fsize = _se[fsize]
local coef_marr = _b[marr]
local se_marr = _se[marr]
local coef_twoearn = _b[twoearn]
local se_twoearn = _se[twoearn]
local coef_db = _b[db]
local se_db = _se[db]
local coef_pira = _b[pira]
local se_pira = _se[pira]
local coef_hown = _b[hown]
local se_hown = _se[hown]
local coef_cons = _b[_cons]
local se_cons = _se[_cons]

putdocx paragraph
putdocx text ("表4 OLS回归完整结果（因变量：净金融资产）"), font("微软雅黑", 10) bold

* 创建三线表
putdocx table tbl4 = (12, 2), border(all, nil) layout(autofitcontents)
putdocx table tbl4(1,.), border(top, single, black, 1.5)
putdocx table tbl4(2,.), border(bottom, single, black, 1)
putdocx table tbl4(12,.), border(bottom, single, black, 1.5)

* 表头
putdocx table tbl4(1,1) = ("变量"), font("微软雅黑", 10) bold halign(left)
putdocx table tbl4(1,2) = ("系数"), font("微软雅黑", 10) bold halign(center)

* 数据行
putdocx table tbl4(2,1) = ("401(k)参与"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(2,2) = ("`coef_e401'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(3,1) = ("年龄"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(3,2) = ("`coef_age'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(4,1) = ("收入（千美元）"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(4,2) = ("`coef_inc'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(5,1) = ("教育年限"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(5,2) = ("`coef_educ'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(6,1) = ("家庭规模"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(6,2) = ("`coef_fsize'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(7,1) = ("已婚"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(7,2) = ("`coef_marr'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(8,1) = ("双收入家庭"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(8,2) = ("`coef_twoearn'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(9,1) = ("固定收益计划"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(9,2) = ("`coef_db'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(10,1) = ("个人退休账户"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(10,2) = ("`coef_pira'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(11,1) = ("住房所有权"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(11,2) = ("`coef_hown'***"), font("微软雅黑", 9) halign(center)

putdocx table tbl4(12,1) = ("常数项"), font("微软雅黑", 9) halign(left)
putdocx table tbl4(12,2) = ("`coef_cons'***"), font("微软雅黑", 9) halign(center)

* 表格注释
putdocx paragraph
putdocx text ("注：括号内为稳健标准误。*** p<0.01, ** p<0.05, * p<0.1。"), font("微软雅黑", 8) italic
putdocx paragraph
putdocx text ("样本量：`ols_n'，R² = `ols_r2'。"), font("微软雅黑", 8) italic

putdocx pagebreak

* 第五章：结论与政策含义
putdocx paragraph, style(Heading1)
putdocx text ("五、结论与政策含义"), font("微软雅黑", 18) bold

putdocx paragraph, style(Heading2) spacing(before, 10)
putdocx text ("5.1 主要发现"), font("微软雅黑", 14) bold

putdocx textblock begin
基于双重机器学习的分析，我们得出以下主要结论：

1. 401(k)退休计划参与显著增加个人净金融资产
2. 控制混杂因素后，参与的因果效应约为`ddml_coef'千美元
3. 这一效应在统计上具有显著性
4. DDML方法有效控制了混杂偏倚
putdocx textblock end

putdocx paragraph, style(Heading2) spacing(before, 15)
putdocx text ("5.2 政策建议"), font("微软雅黑", 14) bold

putdocx textblock begin
基于研究结果，我们提出以下政策建议：

1. 鼓励企业推广401(k)等退休储蓄计划
2. 提供税收优惠激励个人参与退休储蓄
3. 加强金融教育，提高退休储蓄意识
4. 对低收入群体提供额外支持
putdocx textblock end

putdocx paragraph, style(Heading2) spacing(before, 15)
putdocx text ("5.3 研究局限"), font("微软雅黑", 14) bold

putdocx textblock begin
本研究存在以下局限：
1. 观测性研究，无法完全排除未观测混杂
2. 样本代表性可能受限
3. 401(k)参与强度和缴费率信息不足
putdocx textblock end

* 保存Word文档
putdocx save "output/sipp1991/sipp1991_ddml_analysis_report.docx", replace

/*------------------------------------------------------------------------------
第六部分：保存结果
------------------------------------------------------------------------------*/

* 保存分析结果数据
save "output/sipp1991/sipp1991_results.dta", replace

* 导出描述性统计到CSV
summarize $Y $D $X
postfile memhold str20 variable mean sd min max using "output/sipp1991/descriptive_stats.csv", replace
post memhold ("净金融资产") (`r(mean)') (`r(sd)') (`r(min)') (`r(max)')
post memhold ("401(k)参与") (.) (.) (.) (.)
post memhold ("年龄") (.) (.) (.) (.)
post memhold ("收入") (.) (.) (.) (.)
post memhold ("教育") (.) (.) (.) (.)
post memhold ("家庭规模") (.) (.) (.) (.)
post memhold ("婚姻状况") (.) (.) (.) (.)
post memhold ("双收入") (.) (.) (.) (.)
post memhold ("固定收益") (.) (.) (.) (.)
post memhold ("个人退休账户") (.) (.) (.) (.)
post memhold ("住房所有权") (.) (.) (.) (.)
postclose memhold

* 导出DDML结果到CSV
preserve
clear
input str20 method coefficient std_error ci_lower ci_upper
"DDML" `ddml_coef' `ddml_se' `ddml_ci_lower' `ddml_ci_upper'
"OLS" `ols_coef' `ols_se' `ols_ci_lower' `ols_ci_upper'
end
outsheet using "output/sipp1991/causal_effect_estimates.csv", comma replace
restore

log close

* 显示完成信息
display as text "分析完成！"
display as text "输出文件保存在 output/sipp1991/ 目录中"
display as text "包括："
display as text "- Word报告: sipp1991_ddml_analysis_report.docx"
display as text "- 图表文件: figures/ 目录"
display as text "- 数据文件: sipp1991_results.dta"
display as text "- CSV文件: descriptive_stats.csv, causal_effect_estimates.csv"
display as text "- 日志文件: sipp1991_ddml_analysis.log"
